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ABSTRACT 

We performed the first numerical simulations of magnetorotational instability from a sub-magnetar- 
class seed magnetic field in core collapse supernovae. As a result of axisymmetric ideal MHD simula- 
tions, we found that the magnetic field is greatly amplified to magnetar-class strength. In saturation 
phase, a substantial part of the core is dominated by turbulence, and the magnetic field possesses 
dominant large scale components, comparable to the size of the proto-neutron star. A pattern of 
coherent chanel flows, which generally appears during exponential growth phase in previous local 
simulations, is not observed in our global simulations. While the approximate convergence in the 
exponential growth rate is attained by increasing spatial resolution, that of the saturation magnetic 
field is not achieved due to still large numerical diffusion. Although the effect of magnetic field on 
the dynamics is found to be mild, a simulation with a high-enough resolution might result in a larger 
impact. 

Subject headings: supernovae: general — magnetohydrodynamics (MHD) — Instabilities — methods: 
numerical — stars: magnet ars 



1. INTRODUCTION 

The explosion mechanism of core-collapse supernovae 
(CCSNe) is still unresolved despite persistent efforts by 
many researchers over several decades. Recent state-of- 
art simulations show that the neutrino heating mecha- 
nism assisted by hydrodynamical instabilities revives the 
accretion shock, i.e., the explosion results in. However, 
estimated explosion energies are smaller than canoni- 
cal value of orde r 10^^ erg (e.g.. iMarek fc Jankal 120091 
iSuwa et al.1l20l0[ ). 

Meanwhile, effects of magnetic field on the explosion 
dynamics has been studied well for the decade. Nu- 
merical simulations assuming a strong poloidal mag- 
netic field (typically 10^^-10^^ G at the pre-collapse 
phase), and rapid rotation in most cases, show that 
magn etic force assists in driving the energetic explosion 
(e.g., Yamada & Sawai"2004'; 'Obergauhnger et al." "20061 
[Burrows et al. 2007; Takiwaki et al. 2009; Sawai et al] 
[2013). The magnetic fields assumed in these simulations 
are so strong that the conservation of magnetic flux dur- 
ing collapse results in the field strength of > 10^^ G for 
the proto-neutron star surface. This is comparable to 
inferr ed surface magnetic field s of magnetar candidates 
(see Woods fc Thom psonll2006 l, for review of magnetars). 

At present, however, the strength of the magnetic 
field at the pre-collapse stage and the origin of strong 
magnetic field of magnetar are very un c ertain . Stel- 
lar evolution simulations by iHeger et al.l ()2005l ). which 
implement Tayler-Spruit dynamo model, show that the 
pre-collapse strength of the poloidal magnetic field in 
15 Mq star is only 10^ G. Meanwhile, recent observa- 
tions report that some OB stars in the main sequence 
stage possess an ~ 1 kG surface magnetic field (e.g., 
I Wade et al.l |2012[). which corresponds to m agnetar-class 
magnetic flux. iFerrario fc Wickramasingh e (2006) car- 



ried out a population synthesis calculation from main se- 
quence stars to neutron stars, assuming the magnetic flux 
is conserved in the post-main-sequence evolution (fossil 
field hypothesis). Their result implies that ~ 10% of 
OB stars have magnetar-class magnetic flux, while the 
majority have ~l-2 orders of magnitude weaker one. 

Even when the pre-collapse magnetic flux is weak, cor- 
responding to the magnetic field of < 10^^ G for the 
proto-neutron star surface, MHD instabilities may am- 
plify it to magnetar-class strength. So far, there have 
been a small number of wo rks focusing on this issue. 
[Thompson fc DuncanI (|1993f ) argued that convective dy- 
namo in proto-neutron stars generates a n iagnetar-class 
magnetic field. Simulations performed by lEndeve et al.l 
(2012) shows that standing accretion shock instability 
amplifies the magnetic field around the proto-neutron 
star surface from ~ 10^^ G to ~ 10^^ G. For rapidly ro- 
tating progenitors, mag netorotational instability (MRI, 
[Balbus fc HI^ 1991) ma y gre atly amplify the mag- 
netic field (Akivama e t al.l [2003f ). Local box simula- 
tions characterizing a post-bounce core show that an 
initially weak magnetic field, ~ 10^^-10^^ G inside the 
proto-neutron star, exponentially grows due to MRI 
(|Obergauling e r et alj [200l iMasada et al.ll2012[ )PI 

In simulations of MRI from a sub-magnetar-class seed 
magnetic field, < 10^^ G, a quite fine grid size compared 
to the scale of the iron core, e.g., ~ 1000 km, is required: 
The wavelength of maximum growing mode around the 
surface of the proto-neutron star is 
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^ In the simulations performed by Obergau linger et al.l (|2009l ). 
although the size of the simulation box is small, ~ 1 km, compared 
to the iron core, the global gradients of physical quantities are taken 
into account, by which they refer these simulations semi-global. 
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where vp^ is Alfven velocity. In order to attain such a 
high resoluti on, local simulations are adopted in previou s 
works (O bergaulinger et all I2QQ9I : iMasada et all I2Q12D . 
However, in local simulations, it is problematic to pre- 
pare a suitable background state: Although an initially 
stationary background is usually used in local simula- 
tions, post-bounce cores are dynamical. Additionally, 
local simulations are incapable of studying the global dy- 
namics. To relieve these problems, global simulations are 
desirable. 

In this letter, we report the first global simulations 
of MRI in CCSNe from a seed magnetic field of sub- 
magnetar-class fiux. To ease high computational cost 
demanded in global simulations, we carried out compu- 
tations in axisymmetry and in limited radial range, 50- 
500 km. Particular attention is paid to, (i) whether MRI 
amplify the magnetic field to magnetar-class strength 
even in the dynamical background of CCSNe, and (ii) 
whether the amplified magnetic field is strong enough to 
affect the dynamics. 

2. NUMERICAL METHODS 

In the simulations, we solve the following ideal MHD 
equations, employin g a time-exp licit Eulerian MHD 
code, Yamazakura Sawai et al.ll2Q13>) : 
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in which notations of the physical variables follow cus- 
tom. Here, ^ is Newtonian mono-pole gravitational po- 
tential. Computations are done with polar coordinates 
in two dimension, assuming axisymmetry and equato- 
rial symmetry. We a dopt a tabu l ated nuclear equation 
of state produced by iShen et all ()1998 a.b). Neutrinos 
are not dealt with in our simulations. The electron frac- 
tion, which is required to obtain the pressure of a fiuid 
element from the EOS tabl e, is given by a prescription 
suggested by Lie bendorfer^ (l2 005l ). We should note that 
the Liebendorfer's prescription is only valid until bounce. 
iSawai etld] ()2013f ) discussed that adopting this prescrip- 
tion in the post-bounce phase may lead to underestima- 
tion of pressure by upto 20%. 

Before simulating MRI, we first follow the collapse of 
a 15 Mq progenitor star (S. E. Woosley 1995, private 
communication) for the central region of 4000 km ra- 
dius, until ~ 100 ms after bounce (basic run). After the 
central density reaches 10^^ g cm~^, the core is covered 
with Nr X No = 720 x 30 numerical grids, where the spa- 



tial resolution is 0.4-23 km. The core is assumed to be 
rapidly rotating with the initial angular velocity profile 
of 



n{r) = Qo- 



(6) 



where ro = 1000 km and l^o = 3.9 rad s~^, corresponding 
to a millisecond proto-neutron star after collapse. The 
same dip ole-like magnetic fi eld configuration as one em- 
ployed bv I Sawai et al.l ()2013f ) is initially assumed with the 
typical field strength of ~ 10^^ G around the pole inside 
a radius of 1000 km, with which the strength of ~ 10^^ G 
is obtained for the proto-neutron star surface. The ini- 
tial rotational energy and magnetic energy divided by the 
gravitational binding energy are 5.0x 10~^ and 5.3x 10~^, 
respectively. 

To conduct high resolution global simulations to cap- 
ture the growth of MRI in the core, we restrict the nu- 
merical domain to a radial range of 50 < (r/km) < 500 
(MRI run). The initial condition of MRI run is set by 
mapping the data of basic run onto the above domain, 
at 5 ms after bounce. At that time, the density at the 
inner- most grids is (1.5-2.2) x lO^^g cm~^. The inner and 
outer radial boundary conditions for MRI run are given 
by the data of basic run, except for Br. The boundary 
values of Br are given so that the divergence free condi- 
tion of the magnetic field is satisfied. The grid spacing 
is such that the radial and angular grid sizes are same, 
viz. Ar = rA^, at the inner and outer most cells. We 
perform MRI runs with four different grid resolutions, 
where the grid sizes of the inner most cells, Armin, (and 
the numbers of grids, Nr x Nq,) are 12.5 m (8900 x 6400), 
25 m (4700 x 3200), 50 m (2300 x 1600), and 100 m 
(1200 X 800). Note that with above parameters, typi- 
cal MRI growth wavelength around the inner boundary 
is several 100 m at the beginning of MRI runs. We stop 
each MRI run, before the shock surface reaches the outer 
boundary. 

3. RESULTS 

In each MRI run, we found a clear exponential growth 
of the poloidal magnetic energy during t ^ 4-18 ms, 
where t = ms corresponds to the beginning of MRI 
run (see left panel of Figure [1]). The growth timescale is 
larger for higher grid resolution, but almost converges to 
TMRi,num ^ 8 ms at MRI run with Armin = 25 m. After- 
ward, the poloidal magnetic energy gradually increases, 
and roughly saturates around t ^ 50 ms, in which the 
saturation energy is larger for higher grid resolution. 
The evolution of toroidal magnetic energy does not show 
clear exponential growth in each MRI run. Although the 
growth rate is initially larger compared with that of the 
poloidal magnetic energy, it becomes smaller after the 
exponential growth of the poloidal magnetic energy sets 
in. At saturation, in MRI run with Armin = 12.5 m, the 
poloidal magnetic energy is twice as large as the toroidal 
magnetic energy. 

We also follow the variation of the poloidal magnetic 
field structure, for MRI run with Armin = 12.5 m. 
Around t ^ 4 ms, when the exponential growth begins, 
magnetic field lines around the pole in the vicinity of 
the inner boundary start to bend. During the exponen- 
tial growth phase, the magnetic field lines in this region 
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Figure 1. Evolutions of the poloidal magnetic energies (solid 
lines) and the toroidal magnetic energies (dotted lines) for MRI 
run with Armin = 12.5 m (red lines), Armin = 25 m (green lines), 
Armin = 50 m (blue lines), Armin = 100 m (magenta lines), and 
basic run (black lines). 

are further stretched, and filaments of strong magnetic 
field, ~ 10^^ G, appear (see panel (a) of Figure [2] for 
t = 11 ms). After the exponential growth ceases, the 
topology of magnetic field lines becomes rather tangled, 
which implies that the fiow becomes turbulent (see panel 
(b) for t = 20 ms). Subsequently, filaments of strong 
magnetic field and the turbulent fiow pattern, which first 
appear around the pole, also emerge in other regions. 
This phase corresponds to the duration of the gradual 
increase of the poloidal magnetic energy, t ^ 18-50 ms, 
observed in Figure [TJ After the poloidal magnetic energy 
saturates {t ~ 50 ms), a considerable fraction inside a ra- 
dius of 150 km is dominated by turbulent fiow, where the 
strength of the magnetic field reaches ^ 10^^-10^^ G in 
filamentous fiux tubes (see panel (c) for t = 70 ms). 

The bent magnetic field lines, the exponential growth 
of the poloidal magnetic field, the fact the growth rate is 
larger for higher gird resolution, and the turbulence after 
the exponential growth, all invoke the occurrence of MRI 
in our simulations. As seen in panel (e) of Figure [21 the 
number of grids covering the maximum growing wave- 
length in MRI run with Armin = 12.5 m is more than 20 
in most areas, although it is less than five in some limited 
locations in the vicinity of the pole and equator. Thus, 
our highest-resolution MRI run seems almost capable to 
fairly capture t he lin ear growth of MRI. According to 
lAkivama et al.l (|2003l ), the growth timescale obtained by 
the linear theory is 
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In panel (f) of Figure O the distribution of rMRi,th for 
MRI run with Armin = 12.5 m at 4 ms is depicted. It is 
found that, except for the very vicinity of the pole, where 
the grid resolution is not enough, the growth timescale is 
shortest, rMRi,th ~10 ms, around the pole in the vicinity 



of the inner boundary. This is consistent with the fact 
that magnetic field amplification first occurs there. Ad- 
ditionally, the numerically estimated growth timescale, 
TMRi,num ^ 8 ms (Figure [1]), approximately coincides 
with the above theoretical value. We also estimated 
TMRi,th without the effect of the gradient of entropy and 
the lepton fraction, by setting A/" = in Equation ([7j), 
but found insignificant difference from the original one, 
which indicates that convection does not much contribute 
to amplify the magnetic field. With all the above facts, 
we conclude that amplification of poloidal magnetic field 
found in our simulations is mainly due to MRI. 

Since axisymmetry is assumed in our simulations, the 
toroidal component of the magnetic field is not directly 
amplified by MRI. The growth of the toroidal magnetic 
energy, which is slower than that of the poloidal magnetic 
energy, seems due to winding of magnetic field lines by 
differential rotation. 

In lo cal simulations performed bv lObergaulinger et al.l 
()2009f ). the pattern of coherent channel fiows appears 
during the exponential growth phase. In our simulations, 
on the other hand, such a pattern is never observed. Al- 
though chanel-like configurations are locally observed in 
the exponential growth phase (see panel (d) of Figure [2]) , 
they are incoherent, and soon disrupted in timescale of 
milliseconds. Such immediate disruption of channel fiows 
seems to be caused by dynamical fiow motions in the 
background, which is difficult to be taken into account 
in local simulations. 

In left panel of Figure O the energy spectra of the 
poloidal magnetic energy are plotted for MRI run with 
Armin = 12.5 m at different evolutionary phases. Each 
spectrum is derived for radial wave numbers, /c^, and for 
a region 50 < (r/km) < 100 and 15° < 6> < 60°. At 
t = 4 ms, the spectrum shows the dominance of a large 
scale, > 50 km, just reflecting the structure of the back- 
ground magnetic fleld. During the linear growth phase, 
t ~ 4-18 ms, smaller scale structures, < 10 km, grow fast 
compared with the larger ones. Consistently, we found 
that the maximum growing wavelength of MRI around 
the pole is generally ~ 0.1-10 km during this phase. At 
the end of the simulation (saturation phase), the spec- 
trum shows almost flat distribution for r > 5 km, and 
steep decay for r <1 km. In the range between them, a 
slope close to ex k~^^^ is observed, which seems to cor- 
respond to the inertial region of Kologomorov's theory. 
Such spectrum affirms that the ffow is turbulent. Also, 
it indicates that large scale components of the magnetic 
fleld, comparable to the size of the proto-neutron star, 
are produced in our simulations. Since the maximum 
growing wavelength of MRI is generally smaller than this 
scale, as mentioned above, the large scale components are 
likely to be a result of inverse cascade. 

We also compare the average spectra of the poloidal 
magnetic energy for t = 65-70 ms among MRI runs 
(right panel of Fig [3|) . It is found that in a lower res- 
olution run, the spectrum tends to turn the steep decay 
at larger scale, which seems to be due to a larger numeri- 
cal diffusion. It is likely that this causes smaller spectral 
energy for a low resolution run over whole scales, and 
smaller saturation level observed in Figure [H 

Figure m shows that magnetic pressure reaches > 10% 
of matter pressure in some locations, which implies that 
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Figure 2. (a)-(c): Distribution of poloidal magnetic field strength at (a) t = 11 ms, (b) t = 20 ms, and (c) t = 70 ms. Only a part of the 
numerical domain is depicted in each panel, (d): Distribution of poloidal magnetic field strength (color) together with velocity directions 
(arrows) at 11 ms for 50 < (r/km) < 70 and 10° ^ < 30°. (e): Distribution of the maximum growing wavelength of MRI divided by 
the grid scale at t =4 ms. (f): Distribution of the maximum growing timescale of MRI (Equation [7|) at t =4 ms. In panel (e) and (f), 
white-colored areas include the location stable to MRI. All panels are depicted for MRI run with Armin = 12.5 m. 
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Figure 3. Left: Spectra of the poloidal magnetic energy for MRI run with Armin = 12.5 m at t = 4 ms (red hne), t = 11 ms (green hne), 
t = 20 ms (blue hne), and t = 70 ms (magenta hne). Right: Average spectra of the poloidal magnetic energy over t = 65-70 ms, for MRI 
run with Armin = 12.5 m (red line), Armin = 25 m (green line), Armin = 50 m (blue line), and Armin = 100 m (magenta line). Each 
spectrum is derived for radial wave numbers, kr, and for a region 50 < (r/km) < 100 and 15° < 6 < 60°. 
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Figure 4. Distribution of the ratio of magnetic pressure to matter pressure for MRI run with Armin = 12.5 m at t = 70 ms. 
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the magnetic field mildly affects the dynamics. Note, 
however, that we did not obtain the convergence of the 
saturation magnetic field, and a larger saturation level 
is attained for a higher resolution run (Figure [T]). This 
suggest that the simulation with a high-enough resolu- 
tion might result in a larger impact of the magnetic field 
on the dynamics. 

4. DISCUSSION AND CONCLUSION 

We have carried out the first global simulations of MRI 
in CCSNe from a sub-magnetar-class seed magnetic field, 
~ 10^^ G around the surface of the proto-neutron star, 
and rapid rotation, in which two-dimensional ideal MHD 
equations are solved with assumptions of axisymmetry 
and equatorial symmetry. As a result of computations, 
we found that MRI greatly amplifies the magnetic field 
to magnetar-class strength, ~ 10^^-10^^ G, and that the 
magnetic field at saturation phase is dominated by large 
scale structure. Although, the impact of the magnetic 
field on the dynamics found to be mild, we do not obtain 
the convergence of the saturation magnetic field strength, 
which is larger for a higher resolution run. A simulation 
with a higher resolution is necessary to assess the actual 
impact. We also found that the evolution of fiow pat- 
tern in our global simulations are quite different from 
those appear in local simulations bv lObergaulinger et al.l 
(120091) . 

Since MRI and turbulence are intrinsically non- 
axisymmetric phenomena, three-dimensional simulations 
are, in fact, necessary to lead to conclusive results. 
Note, however, that even the world's best comput- 
ers may only capable to simulate global MRI from a 
magnetar-class se ed magnetic field in three dimension. 
lOber gaulinger et al. (2009) found that the main differ- 
ence between two-dimensional and three dimensional 
simulations is that the former results in a larger satu- 
ration magnetic field compared with the latter. They 
argued that the smaller saturation level in three dimen- 
sion is caused by the disruption of coherent chanel fiows 
before they grow well due to non-axisymmetric parasitic 
instabilities. Such difference might not appear in global 
simulations, since our simulations show that coherent 
chanel fiows do not develop even in two dimension. 

While the present simulations set the position of the 
inner boundary at r = 50 km, we found in basic run 
that most rotational energy is reserved inside a radius of 
50 km. It is expected that a simulation with the inner 
boundary located at a smaller radius would result in a 
larger saturation magnetic field. What we obtained in 
this study may be a lower limit. 

Although, in our simulations, a sub-magnetar-class 
magnetic field is amplified to magnetar-class strength 
due to MRI, we should be cautious in concluding that 
MRI can be the origin of magnetar fields. Before that, 
further investigations may be required on the sustain- 
ability of a large scale, strong magnetic field until, e.g.. 



the formation of the neutron star. Additionally, it is 
worth investigating whether magnetar-class strength is 
also attained by MRI from a weaker seed magnetic field, 
e.g., < 10^^ G around the proto-neutron star surface. 
Note that previous local simulations of MRI, includ- 
ing ones in the context of accretion disks, show that a 
weaker see d magnetic field results in a lower saturation 
level (e.g. jOber gauling er et al.ll2QQ9l : lOkuzumi fc Hiro"s3 
I2OIII ). The dependence of the saturation magnetic field 
on the rotational velocity, which also has been discussed 
in many local simulations (e.g., Masada et al. 2012), is 
another issue to be studied in future works. 
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21840050, 24244036.) 
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